Identification of risk areas for Orobanche cumana and Phelipanche aegyptiaca in China, based on the major host plant and CMIP6 climate scenarios

Abstract Parasitic broomrape of the genus Orobanche poses a formidable threat to producing many crops in Europe, Africa, and Asia. Orobanche cumana and Phelipanche aegyptiaca are two of China's most destructive root parasitic plants, causing extreme sunflower, tomato, melon, and tobacco damage. However, the potentially suitable areas of O. cumana and P. aegyptiaca in China have not been predicted, and little is known about the important environmental factors that affect their extension. Due to their invasiveness and economic importance, studying how climate change and host plants may affect broomrapes’ distribution is necessary. In the study, we first predicted the potentially suitable areas of the invasive weeds (O. cumana and P. aegyptiaca) and their susceptible host plants (Helianthus annuus and Solanum lycopersicon) using MaxEnt. Then, the risk zones and distribution shifts of two broomrapes under different climate conditions were identified by incorporating the distribution of their susceptible host plants. The results highlighted that the potential middle‐ and high‐risk zones for O. cumana and P. aegyptiaca amounted to 197.88 × 104 km2 and 12.90 × 104 km2, respectively. Notably, Xinjiang and Inner Mongolia were the highest‐risk areas within the distribution and establishment of O. cumana and P. aegyptiaca. Elevation and topsoil pH were the decisive factors for shaping O. cumana distribution; precipitation seasonality and annual precipitation were the dominant bioclimatic variables limiting the spread of P. aegyptiaca. The potentially suitable areas and risk zones of O. cumana would decrease significantly, and those of P. aegyptiaca would fluctuate slightly under future climate change scenarios. Overall, our study suggested that the two broomrapes’ risk zones will significantly northward to higher latitudes. The results will provide suggestions for preventing O. cumana and P. aegyptiaca.


| INTRODUC TI ON
Parasitic plants are ubiquitous species in ecosystems, including 292 genera and ca. 4,750 species (Nickrent, 2020). Some parasitic plants are considered pathogens due to their negative impact on host plants, such as Orobanche L., Cuscuta L., and Striga L. (Nickrent & Musselman, 2004). Orobanche L. (broomrape) was the typical obligate holoparasitic plant (Gevezova et al., 2012). They depend entirely on the host for survival, which reduces the host's ability to grow and may eventually cause the host's death (Musselman, 1980). Broomrapes include over 100 species and represent some of the most devastating parasitic plants, posing a significant challenge to agricultural production (Bari et al., 2019;Das et al., 2020;Wu & Qiang, 2006). In 2018, it was estimated that global annual losses owing to broomrapes damage were $1.3-2.6 billion (¥ 8.2-16.5 bn) (Ahmad et al., 2018). Orobanche cumana and Phelipanche aegyptiaca are widely distributed in northern China and cause catastrophic damage to tomatoes, sunflower, tobacco, legumes, carrot, and other crops (Wickett et al., 2011). The annual occurrence area only in Xinjiang was about 530-670 km 2 , causing economic losses of more than ¥ 0.5 bn (Yao et al., 2017). The damage induced to the crop by broomrape parasitism can range from 0% to 100%, which varies in each broomrape-host combination (Dhanapal et al., 1996).

Orobanche cumana mainly parasitizes sunflower (Helianthus annuus).
It was discovered for the first time parasitizing cultivated sunflower in Heilongjiang in 1959. Since then, O. cumana has spread over the sunflower cultivation area. In 2020, O. cumana occurred widely in the main sunflower producing area in China, and the occurrence rate was as high as 76.9%, resulting in some fields cannot continue to plant sunflowers, which has become the main factor restricting the sustainable development of sunflowers (Wu et al., 2020). Compared with O. cumana, Phelipanche aegyptiaca could parasitize more hosts (Kacan & Tursun, 2012;Parker, 2012). The tomato (Lycopersicon esculentum) is one of the significant vegetable crops in China that is highly vulnerable to infest by P. aegyptiaca (Tokasi et al., 2014).
In 2013, the damaged area of P. aegyptiaca in tomato had reached about 70 km 2 in Xinjiang province (about 10% of the total tomato planting area), and it spreads rapidly at a rate of 10%-20% every year (Chai et al., 2013;Zhang, Gan, et al., 2016). On average, tomato yield losses were between 30% and 80%, some of which were up to100% for being infected by broomrapes (Das et al., 2020).
It is challenging to control broomrapes because of their concealed underground for the major part of their life, lack of photosynthesis, and physical connection with the crops (Foy et al., 1989).
The problem is exacerbated by the seed-setting rate of broomrape is high and adaptable. Furthermore, this weed has strong invasive potential because the minute seeds can spread rapidly to distant fields in a short time. With climate change, there is possible for this weed to extend its distribution further and become an increasingly severe weed in certain areas. However, until now, no effective prevention measures have been developed for this weed (Libiaková et al., 2018).
An alternative method is to carry out strict plant quarantine to prevent further dispersal in countries where O. cumana and P. aegyptiaca have not yet appeared (Hiroaki & Yukihiro, 2018). Therefore, it is vital to predict the potentially suitable areas of broomrapes and identify their critical environmental variables, which can alert the government to the potential risks on agriculture with future climate change and help them to propose positive coping strategies to prevent the spread of these weeds.
Environmental variables are biologically significant for determining the climatic niche of a species (Mohammadi et al., 2019;Qin, Liu, et al., 2017). Future climate change may have seriously affected global ecosystems and biotas by impacting species distributions (Hulme, 2017;Thomas et al., 2004). In the background of climate change, more and more attention has been concentrated on projecting the future distribution of species to mitigate the influences of future climate change on biodiversity (Ma & Sun, 2018). Species distribution models (SDMs) based on niche theory are a classical approach that uses species' occurrence data combined with bioclimatic data to study the influences of climate change on species distribution (Schipper et al., 2014). Typical SDMs include CLIMEX, GARP, BIOCLIM, DOMAIN, and the maximum entropy (MaxEnt) algorithm (Byeon et al., 2018). Among which, the MaxEnt algorithm is better than other methods in determining potential species distributions based on environmental variables, either in terms of the quantity or in terms of the distribution scope of species records (Hernandez et al., 2010;Phillips & Dudík, 2008).
At present, there has been much research on forecasting the changing trend of potentially suitable areas in some species, such as invasive creatures, dominant forest species, precious and endangered species, but only a few studies have focused on parasitic plants (Ren et al., 2020;Wang et al., 2019;. Climate change can impact the distribution of parasitic plants by changing the geographical regions of host plants. The host plant distribution in conjunction with climate factors can enhance the prediction accuracy of the model (Berzitis et al., 2014). For example, the potential expansion of the co-occurrence between pests and host species increased the harmful implications of fall armyworm (Zacarias, 2020). More recently, when analyzing factors affecting the distribution of emerald ash borer, it was found that adding host plant distributions in models could help with risk assessment (Dang et al., 2021). In parasitic plants, the availability of their main nutrient sources could be another determinant of broomrapes' distribution under future climate. Thus, future control of the broomrapes should consider the influences of climate change and host plant distribution on broomrapes' potentially suitable areas.  (Sauerborn et al., 2007). Orobanche cumana and P. aegyptiaca are the most problematic and damaging broomrape species in China. Botanically, O. cumana is characterized by unbranched inflorescence stems, while P. aegyptiaca is branched (Joel, 2009). The damage of O. cumana in sunflowers has lasted for nearly 100 years (Parker, 2009 Orobanche cumana had been a severe problem to limit sunflower productivity (Shi et al., 2015). Phelipanche aegyptiaca is a highly troublesome parasitic weed in the tomato field (Hershenhorn et al., 2009), and which was first recorded in Aksu of Xinjiang in 1953 (Zhang & Jiang, 1994). Since then, P. aegyptiaca has become the major problem in tomatoes and can totally destroy tomatoes under heavy infestation in Xinjiang.

| Occurrence records collection and screening
The current distribution records for two weeds and their host plant were obtained from among the following aspects: (1)  lycopersicum were obtained for the modeling. All current data were sorted by species name, longitude, and latitude and stored in a CSV format ( Figure 1).

| Environmental variables
A total of 34 environmental variables were screened that may impact the distribution of broomrapes and host plants (  (Fischer et al., 2008). One topographic variable represented by elevation was downloaded from the World Soil Database (http://www.fao.org/faost at/en/#data.) (New et al., 1920).
Future climate scenarios were generated based on downscaled global climate model (GCM) data from CMIP6 (IPPC Sixth Assessment).
Due to the climatic stronger warming, the CMIP6 model differs from CMIP5; the CMIP6 model has higher climate sensitivity and has new specifications in concentration, emission, and socio-economic development (Ramasamy et al., 2021). In the sixth IPCC report, a new set of scenarios-Shared Socioeconomic Pathways (SSPs) have been chosen to drive climate models for CMIP6. In precise, a set of scenarios were selected to support a series of different end-ofcentury climate change results. These scenarios were called SSP126, SSP245, SSP460, and SSP585, each of which represents the distinct probable future greenhouse gas concentration, thus making this set of scenarios ideal for exploring various warming pathways (Gidden et al., 2018). The extreme weather variations may greatly influence the future distribution pattern and its temporal and spatial variations  (Wickham, 2009;Zhao et al., 2021). The environmental variables with absolute values of the Spearman correlation coefficient of less than 0.8 correlated ecological physiology and could be used for modeling (Qin et al., 2019). Then, the VIF analysis was performed on all variables, and the environmental factors with a VIF <10 were selected (Ye et al., 2020)

| MaxEnt modeling approach, optimization, and evaluation
In this research, the potential geographical distribution of broomrapes (O. cumana and P. aegyptiaca) and host plants (H. annuus and S. lycopersicum) under current and various future climate conditions were carried out in MaxEnt (maxent 3.4.1; http://www. cs.princ eton.edu/wscha pire/maxent) with presence-only data.
For each species, the training data were 75% of the sample distribution points chosen randomly, and the testing data were the remaining 25% of the sample points. The MaxEnt was run with a maximum number of iterations of 10,000, the 10 −5 default convergence threshold, and 10 replicates under bootstrap run type to guarantee the model's accuracy . To ensure that the possibility of these species' distribution appeared close to normal, the package "ENMeval" was used in R x64 4.0.2 (http:// Jackknife testing was applied to assess the importance of variables. The relative importance of a single factor was evaluated by using default estimates of percentage contributions (Elith et al., 2015). Receiver operating characteristic (ROC) analysis was applied to assess the performance of the model, and the area under the ROC curve (AUC) was applied as an excellent index to provide the accuracy estimates (Katz & Zellmer, 2018;Peterson et al., 2017;Vanagas, 2004). Usually, AUC values were between 0 and 1. AUC < 0.5 represents a model that predicts worse than random and rarely occurs in reality, while AUC > 0.75 describes a model with a fair ability to discriminate (Ferrier, 2000). The closer the AUC is to 1, the better the prediction ability is (Yi et al., 2016). According to the AUC, prediction ability is classified as failing (0.5-0.6), poor (0.6-0.7), fair (0.7-0.8), good (0.8-0.9), or excellent (0.9-1) (Swets, 1988).

| Analysis of spatial pattern changes
The spatial units with a species presence possibility value ≥0.3 was regarded as a suitable area for the species, and the value <0.3 was considered an unsuitable area. According to the (0,1) matrix table, the changes in the spatial pattern of the potential risk zones of two weeds under the different climate conditions were further analyzed.
The area changes in the future were computed based on the present risk zones of the species. Four types of changes in the risk zones for species were defined: matrix values 0 → 1 were newly expanded suitable areas, 1 → 0 were lost suitable areas, 1 → 1 were reserved suitable areas, and 0 → 0 were unsuitable areas (Ye et al., 2020).
Lastly, the matrix change value was loaded into ArcGIS 10.4.1 to realize the visualization of the spatial pattern change in different periods.

| Centroids shifts
The changes in risk zones were computed, and the core distributional were compared for current and future risk zones by utilizing the SDM Toolbox (an ArcGIS tool) in ArcGIS10.4.1 (Brown, 2014).
The focus of this analysis was to summarize the core distribution shifts in the risk zones of these two weeds. This was achieved by reducing the species' risk zones to a center point and establishing a vector particle that reflected the magnitude and direction of the predicted changes through time. Finally, the shifts in species' risk zones were checked by tracking how the core distributional changed with different SDMs, and the "dists" function in R x64 4.0.2 was used to calculate the centroid migration distance of the risk zones in latitude and longitude coordinates (Zhang et al., 2018).

| Species distribution model performance
According to occurrences record of species and environmental data, we obtained potential distribution maps where broomrapes and host plants might occur. Models for the broomrapes and host plants performed better than random with the given series of training and test data. The average values of training AUC and testing AUC for all models were greater than 0.8, suggesting that these models were accurate for all studied species predictions and could be reliably used to simulate the relationship between diversity patterns and climatic variables ( Figure S2).
Using the response curve ( Figure S3), we gained the thresholds

| Potential geographical distribution under current climate conditions
According to occurrence records and environmental data, the current potentially suitable areas of two broomrapes are illustrated in   (Figure 2c).
The potentially suitable areas of H. annuus and S. lycopersicon based on current environmental factors are shown in Figure   S4 and According to the prediction results from MaxEnt, the potentially suitable areas of broomrape and its host plant were superimposed to obtain the risk zones of broomrape (Figure 2b 197.88 × 10 4 km 2 , is currently a potential middle-and high-risk zone (Table 2). For P. aegyptiaca, the middle-and high-risk zones were detected in Xinjiang (Figure 2d). The predicted area was around 12.91 × 10 4 km 2 , 1.34% of China's area (Table 2).

O. cumana
The calculated results for the O. cumana habitat suitability in SSP126, SSP245, and SSP585 are illustrated in Figure 3-Ⅰ, Ⅲ , and

P. aegyptiaca
The distribution maps of suitable areas for P. aegyptiaca in future climate conditions for 2050 and 2090 are summarized in Figure 4-Ⅰ, Ⅲ , and Table 2. Compared with current climate conditions, the total suitable area of P. aegyptiaca in SSP126 and SSP585 was projected to decline slightly but increase slightly under the SSP245 scenario ( Figure S6-b). Highly suitable areas were still concentrated in northern Xinjiang, and the suitable areas ranged from 20.52 × 10 4 km 2 -39.69 × 10 4 km 2 , occupying 2.13%-4.12% of China's area.
The predicted results of future areas' suitability under different SSPs are shown that the changing climatic condition alters the distribution of suitable future areas of S. lycopersicum ( Figure S5-Ⅲ, Ⅵ and Table S3). Compared to the current potential distribution, for 2050, the slight increase in area for high suitability was found under SSP126, and a slight decrease was found under SSP245 and SSP585.
In 2090, the highly suitable areas of S. lycopersicum increased with increasing temperature. Under SSP585, the high suitability area increased significantly to 193.75 × 10 4 km 2 .   Under the SSP126-2050 and SSP585-2050 scenarios, MaxEnt predicted that P. aegyptiaca total losses would occur mainly in western Xinjiang, with 0.07 × 10 4 km 2 and 3.36 × 10 4 km 2 , respectively.

| Future changes in risk zones
Conversely, in SSP245-2050, the total risk habitats expanded by 4.52 × 10 4 km 2 , and the new areas were widely concentrated in northern Xinjiang, which indicated that a lot of the low-risk areas under this condition might be replaced by the middle-risk suitable areas ( Figure 5-Ⅲ and Table S4). By 2090, the MaxEnt model predicted total losses of 6.51% and 12.17% of the current area under the SSP126 and SSP245 scenarios, respectively. The distribution of newly added risk zones was also concentrated in northern Xinjiang under the SSP585-2090s combination but was slightly less than the scenario of SSP245-2050s at an increase of 4.37 × 10 4 km 2 (Figure 5-Ⅳ and Table S4).  For the risk zones of Beijing, the core of the current risk zone was located at 115°47′E and 39°59′N (Figure 6-b). Under all future climate scenarios (SSP126, SSP245, and SSP585), we saw that the cen-

| DISCUSS ION
Broomrapes have been listed as entry plant quarantine pests and agricultural plant quarantine pests in China (Wu & Qiang, 2006). temperatures required for its conditioning and germination, which was found to be around 25-30°C: both lower and higher temperatures resulted in poor germination (Shi et al., 2018). Meanwhile, this also may be attributed to transitory or constantly high atmospheric temperatures causing a series of heat injuries, physiological disorders, and morphological and biochemical changes in plants (Lipiec et al., 2013).
Climate is a critical factor in predicting the potential distribution of species. However, the distribution of susceptible host plants is  (Johkan et al., 2011). This might reduce the risk zones of O. cumana.
In addition, previous studies suggested that temperature impacts on the host-parasite relationship are complex (Sukno et al., 2001).
Sukno et al. studied the effects of temperature on the resistance of sunflower to parasitism by O. cumana and found that the level of parasitism was lowest for all the treatments at the highest temperature (Sukno et al., 2001). These results have certain similarities with those of  and , who found that under controlled conditions in polyethylene bags and Petri dishes, as the temperature rose, more Orobanche tubercles degenerated and died; Namely, the sunflower plants expressed higher levels of resistance. The defense mechanism may have increased with rising temperature, causing the broomrapes to develop retardation, necrosis, or die. Taken together, the suitability of conditions for the host can also influence the future suitability of conditions for parasites.
However, incorporating biotic interactions (predation, competition, parasitism, and mutualism) in SDMs does not always improve predictive accuracy (Pellissier et al., 2010;Silva et al., 2014), as it requires extensive knowledge regarding biotic interactions and populationlevel impacts, as well as detailed and specific information in biotic interactors and target species (Fordham et al., 2012;Simões & Peterson, 2018 (Hu et al., 2015). During the broomrapes' life cycle, the temperature is the dominant factor affecting their growth, and there is a strong interdependence between parasitic ability and temperature (Eizenberg & Goldwasser, 2018). Ephrath and Eizenberg (2010) reported that the relationship between the process of O. cumana and P. aegyptiaca attachment and development and thermal time followed a sigmoid-shaped curve, and the broomrapes' biomass was also found to be closely correlated to thermal time. Research on the influences of temperature on the germination of different species of broomrapes showed that each broomrape had a specific best temperature range for germination and growth, which generally reflected its geographical distribution (Song et al., 2005). The optimal temperatures for both conditioning and germination were about 25-30°C for O. cumana and about 18-21°C for P. aegyptiaca (Shi et al., 2018;Song et al., 2010). Final germination must fulfill the following two conditions: the minimum average temperature that must be exceeded and the maximum temperature above which the seeds will not germinate (Ermias & Murdoch, 1999). Therefore, appropriate adjustment of the sowing date can considerably reduce the abundance of broomrapes and lessen the damage caused by broomrapes to host plants.
Humidity is another major factor affecting broomrapes seed germination. The germination of broomrapes seeds requires preculturing under specific temperature and humidity conditions (Kebreab & Murdoch, 1999). When the soil humidity is too high, the parasitism and growth of broomrapes are inhibited, which may be because the humidity affects the dormancy characteristics of broomrapes, and then affects their germination and parasitism (Kebreab & Murdoch, 2000). Less infestation due to broomrapes has been observed in sunflower/tomato grown under flooded irrigation or in years with more precipitation (Punia, 2014). In addition, a substantial reduction in seed viability/longevity was observed when the seeds were imbibed and at higher temperatures. Thus, combining irrigation with solarization during the dry season could significantly raise the depletion of the soil seed bank.
Soil conditions also have a significant influence on shaping the distribution range of broomrapes. Alkaline soil (pH > 7.0) was found to be favorable for broomrape survival and growth, but the parasitism rate was less in acidic soil (pH < 7.0) (Dinesha & Dhanapal, 2013;Di et al., 2020). Shi et al. (2018) reported that sandy loamy soil with a temperature of around 25-30°C, moisture of 60%-70%, and soil pH value of 8 could benefit the attachment and growth broomrapes. Some research results have revealed that the parasitism rates of broomrapes were closely related to the amount of available soil nitrogen and the type of nitrogen fertilizer used (Irmaileh, 1994;Mesbah et al., 2012). As soil nitrogen rates increased, the numbers and dry weights of broomrape shoots declined, the yields of sunflower/tomato substantially rose, and urea fertilizers could significantly reduce the occurrence of broomrapes (Mariam & Suvvanketnikom, 2004;Narges & Syyed, 2014). Therefore, N should be applied at least at the recommended rate if using N at higher doses is not possible on the ground of affordability by farmers and cost economics. In this study, the CMIP6 climatic scenarios were employed to determine the potential effects of climate change on the environmental suitability for broomrapes. Among the 19 environmental variables used for the developed model, temperature, precipitation, and soil properties all made important contributions to the habitat suitability of O. cumana and P. aegyptiaca, suggesting that these factors play essential roles in broomrapes distribution.
Responses to climate change in biological communities are usually related to species distribution about latitude or extreme elevations (Lenoir et al., 2008). Generally, with the climate changes, suitable areas in the Northern Hemisphere are expected to expand northward, and these areas will become more climatically suitable . Some studies have revealed that some species may expand their ranges to higher latitudes and elevations as global warming continues (Bertrand et al., 2011;Walther et al., 2005;Wilson et al., 2007). The study suggested that a northward shift of broomrapes' risk zones to higher latitudes would gradually become more significant. The centroid of risk zones for O. cumuan and P. aegyptiaca shifted to areas at higher latitudes than their current habitat, which may be facilitated through adaptations. Therefore, the long-term risks are greater in higher latitudes, as these areas will experience warming, which will promote the invasion and establishment of O. cumuan and P. aegyptiaca.

| CON CLUS IONS
In this research, we established a predictive model that considered climate, soil, and topography and provided a basis for assessing the suitable area of Orobanche cumana and Phelipanche aegyptiaca. This model was applied to examine the influences of climate change on the potentially suitable area for those two broomrapes in China. We also took global climate change and host plants into account to examine the potential risk status of the broomrapes in the future. The results showed that Xinjiang and Inner Mongolia would experience the greatest risk of establishing the O. cumana and P. aegyptiaca, which may affect local agricultural productivity or pay a greater price to control these weeds. Although the results suggested that the distribution of risk zones for O. cumuan will be significantly decreased with global warming, prevention and treatment should not be taken lightly due to the rapid spread of broomrapes and the difficulty of eradication. The predicted results can be applied to clearly understand the distribution of suitable areas and risk zones of two broomrapes and determine their expansion reasons to provide a useful reference in formulating quarantine and control measures for O. cumana and P.
aegyptiaca. Future studies could be needed to consider the effects of human behavior or the habitats of other host plants, and even the study area could be extended into the global scope.

ACK N OWLED G M ENTS
This work was financially supported by National Natural Science

DATA AVA I L A B I L I T Y S TAT E M E N T
Climate data, MaxEnt input files and the four species' distribution data: Dryad https://doi.org/10.5061/dryad.8pk0p 2nph; and the data that supports the findings of this study are available in the Appendix S1 of this article.